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*Professor, Old Dominion University. Norfolk, VA 23529-0077, U.S.A. 
**Graduate Student, Old Dominion University, Norfolk, VA 23529-0077, U.S.A. 


Introduction 

Previous accomplishments resulting from the research Grant NCC 1-42 have been 
numerous. These accomplishments are summarized in the appendices A,B and C of this 
report. The Appendix A lists publications that have resulted from the varied research 
efforts associated with this grant. The Appendix B lists graduate students who have been 
associated with this grant, their research efforts and their degree accomplishments. The 
Appendix C is a summary of efforts involving Green’s Functions. 

Current research on Grant NCC 1-42 involves the development of a multigroup method 
for the calculation of low energy evaporation neutron fluences associated with the Boltz- 
mann equation. This research will enable one to predict radiation exposure under a variety 
of circumstances. Knowledge of radiation exposure in a free-space environment is a neces- 
sity for space travel, high altitude space planes and satellite design. This is because certain 
radiation environments can cause damage to biological and electronic systems involving 
both short term and long term effects. 

By having apriori knowledge of the environment one can use piediction techniques 
to estimate radiation damage to such systems. Appropriate shielding can be designed 
to protect both humans and electronic systems that are exposed to a known radiation 
environment. This is the goal of the current research efforts involving the multi-group 

method and the Green’s function approach. 

Reference [1] presents a short history of the study of the propogation of space radia- 
tion through matter, the development of space-radiation physics and protection techniques. 
This reference outlines major radiation studies and their results. The accurate prediction 
techniques for dose fields requires either Monte Carlo type solutions, which require a great 
amount of computational time, or a study of the Boltzmann equation under various cir- 
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cumstances. The Boltzmann equation does not lend itself readily to analytical solutions 
and so various types of numerical solutions have been developed. One such numerical 
solution is the HZETRN code developed by Wilson et.al., reference [1], This code has the 
ability to predict dose fields in simulated tissue behind a shield for high charge and high 
energy particles. Flux predicitions by HZETRN are based upon a straight-ahead transport 
of evaporation neutrons with one dimensional angular transport. In this research we in- 
vestigate the transport of evaporation neutrons through a shield-target environment based 
upon source terms generated by the HZETRN code. 

Formulation of Transport Equations 


We define the differential operator 

B[4>} = 


i~jE s ’ {E)+a = {E) 


<p(x, E) 


8<P ^' E) - JL [ ; Sj{E)4>{x , E)} + aj(E) E) 

ox oE 


(i) 


and consider the Boltzmann equation from reference [1] 

pOC 

B[4>j] = J2 / fjk{EtE f )4> k {x,E')dE' 

u Jo 


(2) 


where <f>j is the differential flux spectrum for the type j particles (j — n neutrons or 
j = p protons), S j (E ) is the stopping power of the type j particles and crj{E) is the total 
cross section. The term fjk{E,E f ) is a macroscopic differential energy cross section for 
redistribution of particle type and energy. We write 


fjk{E, E') = J2PW( E ')fjkAE>E') 

13 

where fj^^E^E') is an elastic collision term, op is a microscopic cross section and pp is 
the number density of (5 type atoms per unit mass. The collision terms are expressed as 


fjk,f3 — !jk,p + fjk,p 

where ^ represents evaporation terms and fj k ^ represents direct cascading terms. The 
evaporation process dominates over the low energies (E < 25 Mev) and the direct cascading 
effect dominates over the high energy range ( E > 25 Mev). The equation (2) is written as 

pOC 

b M = E/ + (3) 

t Jo 3 
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For j = n we write equation (3) in the operator form 


B [4>\ — I e [<t>\ + Id\4>\ + Ie[<fip\ + Idl&p] (^) 

where 0 = 0 n and I e ,I d are integral operators. Let 0d denote the solution to the transport 
equation 

B[<t>d] — Id[4>d\ + h l[4>p\ ^ 

which is valid over all energy ranges. The solution to equation (4) is then 0 = 0 e + 4>d 
where 0 e is the neutron flux due to evaporation and (pd represents neutron flux due to 
direct cascading effects. The equation (4) can be expressed in the form 

B[(p e \ + B[(pd\ = ie[0e] + ^e[<Pd] + J d[4>e ] + Id[<frd\ + Id[4>p\ + [<fip] * (®) 

Using the assumption that the numerical solution to equation (5) is readily obtainable 
from the the HZETRN code, reference [6], and that I d [<f> e ] w 0, the equation (6) reduces 
to 

B[(p e \ — I e [<f>e\ + h[<t>d\ + h[<t>p\ ^ 

The stopping power Sj(E) = 0 for neutrons and so the equation (7) reduces to the integro- 
differential transport equation with source term 

rE/ap 

+ cr(E) <p e (x , E) = y] / / s,p{E i E')cfi e (x , E ) dE +g(E,x) (&) 

J 

which represents the steady state low energy neutron fluence depth x and 

energy E. The various terms in equation (8) are energy E with units of (Mev), depth in 
medium is x with units of (g/cm 2 ), 4> e (x,E) (# particles/cm 2 — Mev) is the fluence and 
g(E t x ) = I e [<t>d\ + h\<t>p] (#particles/g - Mev) is a source term. In addition, the equation 
(8) contains the scattering terms 



fs.p = P30p{E')fj kt p(E, E') 


with units of cm 2 /g - Mev. The limits of integration (E, E/ap) represent cut off values for 
neutron production because secondary neutrons produced have approximately the same 
energy as the projectile primary neutrons. The term ap is defined as the ratio 


a 0 



( 9 ) 
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and is a constant less than 1, where A t { is the atomic weight of the i th type of atom being 
bombarded. The quantity a has units of (cm 2 /g) is a macroscopic cross section given by 

a = ^p j a j ( 10 ) 

j 

where pj is the number of atoms per gram and crj is a microscopic cross section in units 
of cm 2 /atom. The reference [2] provides approximate Maxwellian averages of cross section 
values in barns. These values are listed in Table 1 along with other parameters of interest 
for selected elements. 

Other units for equation (8) are obtained from the previous units by using the scale 
factor representing the density of the the material in units of g/crn 3 . 


Table 1. Parameter Values for Selected Elements 

Element 

Symbol 

At, 

Cross Section 
barns* 

Density 

g/cm 3 

a 

Lithium 

Li 

7 

1.050 

.534 

.563 

Carbon 

C 

12 

4.739 

3.52 

.716 

Aluminium 

A1 

27 

1.348 

2.7 

.862 

Calcium 

Ca 

40 

2.99 

1.54 

.905 

Iron 

Fe 

56 

11.40 

7.85 

.931 

Lead 

Pb 

207 

11.194 

11.342 

.981 

* Maxwellian averages (elastic) 


Mean Value Theorem 

Through out the following discussions we employ the following mean value theorem 
for integrals. 


M eanV alueT heorem For <p(x),f(x) continuous over an interval a < x < b such that 
(i) <j>(x ) does not change sign over the interval (a, 6), (ii) <p(x) is integrable over the interval 
(a, b), and (iii) f(x) is bounded over the interval (a, 6), then there exists at least one point 
£ such that 


f(x)<p{x) dx = /(£) 



a < £ < b 
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Multi-group method 


We consider the case (3=1 which represents neutron penetration into a single element 
and let (p e = 0. We integrate the equation (8) from Ej to £h-bi with respect to the energy 
E to obtain 



E) 

dx 


dE + 



<j(E)(t>{x, E) dE = It + 


(ii) 


where 


h = 




f Si0 {E,E')4>{x,E') dE' dE 


( 12 ) 


and 


ti = 



g{E,x) dE. 


(13) 


As a first approximation test case we use the approximate source and scattering terms 
g = g(E,x) = KEe~ E / T , with K and T constants, and f s ^(E,E') = — L _ e( i — so 
that the equation (13) is easily integrated to obtain 


£ t = KT ( E x e~ BilT - E i+ ie" Bl+l/T ) + KT 2 (e £i/T - e Ei + 1 / T 


(14) 


We define the quantity 


*<(*) 



(15) 


and write the equation (11) in terms of the ^i(x) terms as follows. In the first term in 
equation (11), we interchange the order of integration and differentiation to obtain 



8<p(x , E) 

dx 


dE 


d$i(x) 

dx 


(16) 


Using one of the several mean value theorems for integrals, the second term in equation 
(11) can be expressed as 

fEi+l 

/ o<t>{x,E) dE = a$i(x) (17) 

J Ei 

where o~ = o{Ei + 6(Ei+ \ — E t )), for some value of 6 between 0 and 1. 

For the term Ii in equation (12) we interchange the order of integration as illustrated 
in the figures 1(a) (b). The integration of equation (12) depends upon the energy partition 
selected. For example, the figure 1(b) illustrates an energy partition where 1 < Ei/a 
and in this case we can write the equation (12) as 
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rE i + 1 r E‘ rEJa rE i+l rEi+l/* fEi+l 

h = / / HdEdE'+ / / //d£d£'+ / / H dEdE 

J E' = E{ J E=Ei JE' = E i + x J E=Ei J E' = EJa J E=aE’ 


( 18 ) 


where £ = /,(£, E')^(x, £')• The figure 1(c) depicts the case where £,+1 = Ei/a exactly 
for all i. In this special case the equation (12) reduces to 

rEi^ i rE‘ r&i+ i/ Q fEi~v 


rE' rEi + i/ot rE, + i 

Ii = / / H dEdE' + / H dEdE . 

JE’=E i JE-Ei JE'=E i+ , JE=aE' 


(19) 


The integrand £ can be integrated with respect to E and the results expressed in terms 
of the quantities 


and 


E (b, a) 

G(E') 


/ re rE dE = e Tb — e Ta 

a 

a(E')e- rE ' 


( 20 ) 


1 - e -(l-a )tE‘ 

and we can write equation (19) in the form 

Ii= f G{E')F{E',Ei)<i>{x,E') dE 7 

JE'=-Ei 

+ f E,+l/ G(E , )F(E i j r i,aE')(f>(x, E r ) dE' . 

JE' = E , + 1 

To illustrate the basic idea behind the multigroup method we use one of the many mean 
value theorems for integrals and write the equation (20) in the form 

Ii = G{E*)F(E*, Ei)$i + G(E* +l )F{Ei+i,aE? + Jfci+i- 

where E x < E* < Ei/a and E i+i < E* +l < E i+ i/o. The special partitioning of the energy 
as illustrated in the figure 1(c) enables us to obtain from the equation (11) a system of 
ordinary differential equations 


d 

dx 


■ $0 ■ 


"an 

a 12 





U 22 a 23 


$N-2 



CLN- l,N-l 

aN-l,N 

-$N- 1 - 


_ 


a N N 


-1 

$0 




<*>1 

+ 



$N - 2 


£jV- 2 




-£n-i - 


( 21 ) 


where a iyi = G(E*)F{E *, E { ) -a and aq i+ i = G(E* +l )F{E i+ i , aE* +l ). We further assume 
that for large values of N that = 0 for all i > N. This gives rise to the system of ordinary 
differential equations 

dy r 

— — Ay + b 

dx 
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subject to the initial conditions y(0) = 0. Here y is the column vector 

col($o,$i,...,$JV-i), -4 is an N by N upper triangular matrix and b is the column 
vector col(£o,fi,...,&v-i)- 

In a similar manner the integrals in equation (18) can be evaluated for other kinds of 
energy partitioning and we will obtain a system of equations having the form of equation 
(21). However, for these other energy partitions the structure of the N by N square 
matrix A will change. It remains upper triangular but with more off diagonal elements 
which depend upon the energy partition. For our purposes the system of equations (21) 
will be used to discuss some of the problems associated with the multigroup method. 

We construct the energy partition 

{E 0 , Eo/a, Eo/a 2 , ... , Eo/a N } , 

where Eq = 0.1 Mev, for the selected elements of lithium, aluminum, and lead. The table 
2 illustrates integer values of N necessary to achieve energies greater than 30 Mev. 


Table 2. Energy Partition Size N 

Element 

a 

N 

0.1/a* 

Lithium 

0.563 

10 

31.53 

Aluminium 

0.862 

39 

32.75 

Lead 

0.981 

298 

30.38 


For energy partitions where Ei+i < Ei/ a the values of N will be larger and when 
Ei - |-i > Ei/a. the values of N will be smaller. The cases where E{+ 1 > Ei/a give rise to 
situations like that illustrated in the figure 1(d). In this figure the area A\ is associated 
with the integral defining <lq and the area Aj is a remaining area associated with an 
integral which is some fraction of the integral defining 4>i+i which is outside the range of 
integration and so some approximation must be made to define this fractional part. This 
type of partitioning produces errors, due to any approximations, but it has the advantage 
of greatly reducing the size of the iV by N matrix A. 

The case of neutron penetration into a composite material gives rise to the case where 
(3 > 1 in the equation (8). In this special case the equation (12) becomes 


'«-E 

3 




/^(E^E'^ix^E^dEUE. 
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We select a = max(ai, a- 2 , . . . , aj) and construct the energy partition where Ei + \ = Eija. 
We then obtain a system of differential equations having the upper triangular form 


d 

dx 


1 

1 


: ! 


1 

^ ... . 

1 

l 



an 


a 12 
« 22 


a 13 
a 23 

033 


aiN 

CI-2N 


a-NN 



$0 




$1 




; 

+ 

; 


.^n-i _ 


,£iV-i . 


( 22 ) 


Energy Partition for Finite Values of N 

Consider the case of neutron fluence in a single shield material with the energy parti- 
tioning as illustrated in the figure 1(c). This is the case where successive energy values are 
given by Ei + 1 = E{ja for all values of the index i as i ranges from 0 to N. We select a finite 
value for N , say N = 10, and select Eq large enough such that the assumption = 0 
holds true. The system of equations (21) is then a closed system and we can solve for 
the terms {$ 0 , • • • , ^iv-i}- If we march backwards N energy partitions from the original 
starting value Eq , we obtain a new value for Eq , such that the final value $ n equals the 
old starting value <&o- The term ajv,jv+ i$n in the equation (21) is now known and can 
be moved into the right hand side of the system along with the £n-i term. Continuing 
in this manner we can define groups of energy partitions of size N, where in each group 
Ei+i = Ei/a are the energy values considered. The starting value Eq changes for each 
group and, except for the highest energy group, we will have the value of 4>o from a higher 
group equal to the value of 3>jv from the lower energy group. The grouping scheme is illus- 
trated in the Figure 2. In this way, we can break a large energy partitioning into groups 
of N equations, where the highest energy group of equations is solved first and the lowest 
group of equations is solved last. The nonzero elements a l:J for the matrix A in equation 
(22) consists of the diagonal elements and the first diagonal above the main diagonal. This 
gives the values 

a ii =G(E*)F(E*,E 1 )-a 
ai , i+ 1 =G(E* +l )F(E i+ i,aE* +1 ) 

for i = 1 where E* and E* + 1 are selected mean values associated with the lower 

and upper triangles illustrated in the figure 1(c). These mean values vary with energy and 
were selected as 

E* =E i + 9 l (E i+1 -E i ) 

E * + 1 =E i+l + 0 2 (E i+ 2 - E i+1 ) 
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where 


and 


where 


1 

f 7l + mn(E - En) - 6 1 

E > En 

0 i = \ 

7l + m i 2 (E - En) ~ Si 

E 2 2 < E < En 

1 

{ 73 + mis(E - E 22) - £1 

E < E22 


f 72 + rn 2 i{E - En) 

E > En 

0 2 = 

-- < 72 + m- 2 2 (E - En) 

E 22 < E < En 


{ 74 + m-23(E - E22) 

E < E22 

71 =0.93 mn =0.0030485 

m 2 1 =0.004355 

72 =0.90 m v2 =0.2490258 

m-22 = 0.24 9 0 26 

73 =0.30 mi3 = - 0.3937186 

m 23 = - 0.255920 

74 =0.27 En =3.037829 

E22 =0.5079704 


and <5i has the values, 0.0 for lead, 0.02 for aluminum, and 0.075 for lithium. In addition, 
we set 9 3 = 6 1 and 9 4 = 9*2 because of the symmetry of the triangles involved in the 
calculations. The above values of 9 for the mean value theorems were determined by trial 
and error so that the multigroup curves would have the correct shape. These selections 
for the mean values are not unique. 

Solution Method Shield Material 

We consider the energy partition Ei + 1 — Ei/a and the resulting system of equations 
(21). The solution of this system of equations is obtained by first solving the last equation 
of the system. This equation has the form 

d$N- 1 


dx 


&NN 


$N - 1 +&V-1OO1 $iv-l(0) - 0 


and has the solution 

$*-!(*) = e aNNX (* N -i{0) + jf ^jv-i('S)e -a ' VNS ds 

which implies 

*N-l{xo + Ax) = e aNN *** N -i{x 0 ) + c «™(*o+A*> / ^_i(s)e~ aNN3 ds. 

J XQ 


We then consider each of the remaining equations above the last equation. A typical 
equation from this stack has the form 

d<& i 


dx 


= au$i + fi{x), ^t(0) = 0 
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where /;(x) = £i(x) + aj i , + i.$j^i(x) is known since $; + i(x) is calculated before $j(x). This 
typical equation has the solution 

<*>;(*) = a a “* (^>i( 0) + fi(s)e~ a ' iS ds 


which implies 


$i(xo + Ax) = e a ” Ax <l> t (xo) + e 


,au(:ro+A 


rx c 
s) / 

J XQ 




/i(s)e 


Aji 3 


ds 


Another Viewpoint for the Resulting System of Equations 

The equations resulting from the energy partitioning of each group can be expressed 
as a system of ordinary differential equations 

d<$> 

— — + <f$j = Ii + £i, i = li 2, . . • , N - 1 

ax 


which is equivalent to the vector system of ordinary differential equations 


dx 


subject to the initial condition y{ 0) — 0. Here y is the column vector col($o } • • • > $N- i)» 
A is a N by N upper triangular matrix, and b is the column vector co1(£o 7 £l> • • • 1 ) 

From the solution of this system of ordinary differential equations we calculate the average 
fluence over each energy interval 


^i—avg — 



dx. 


Fundamental matrix solution 


Let Y(x) = e Ax denote the fundamental matrix solution defined by the matrix differ- 
ential equation 

^j- = A Y, Y(0) = I 

dx 

where I is the N by N identity matrix and A is an upper triangular matrix. The solution 
of the system of equations (21) can then be represented in the form 


y(x) = Y{x) 



y _1 (s)6 ds = 



(23) 
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The exponential matrix V'(x) satisfies the properties that 

AY(x) =Y{x)A 

y(-x)=y _1 (x) (24) 

and V'(x + s) —Y (x)T(s) 


Consequently, for b constant, the solution given by equation (24) can be represented in the 
form 

y(x) = {Y{x) - I)A~ l b. (25) 

The solution, as given by equation (25) can also be expressed in terms of a Green s function 
for discrete systems. We rewrite the solution system (25) in the equivalent form 

y(x + h) = {Y (x + h) — I)A l b 
y( x + h)=(Y(x)Y(h)-I)A- 1 b 

( 26 ) 

y(x + h) =Y{h)(Y(x) - Y(-h) + I - I) A l b 
y(x + h) =Y(h)y(x) + {Y(h) - I)A 1 b. 


If b is not constant, but a function of x, then the solution is left in the integral form of 
equation (23) and becomes 


y(x) = V'(x) 





In this case, we have 


f*x+h 


y(x + h) = T(x + h) / Y~ 1 {s)b{s)ds = Y{h)y(x) + Y{x 


rX-hh 

) Y (h) l 


Y l (s)b(s) ds. 


Calculation of the Exponential Matrix Exp(Ax) 

The fundamental matrix solution Y(x) = e Ax can be calculated from the Putzer 
algorithm, reference [3]. This algorithm states that for A an n x n matrix with eigenvalues 
Ai, A- 2 , . . . A n , the exponential matrix is given by 

n — 1 

= ^r i+ 1 (x)Pj 
j=o 


where 

P 0 = I, Pj = n J fc=1 (.4 - Afc/), j = 1, n 
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and ri(x), . . . , r n (x) are solutions of the triangular system 

~7~ = A 1 r l r l(°) = 1 
ax 

p- =rj_d*) + rj(z), rj(0) = 0 
ax 

for j = 2, . . . , n. Here each eigenvalue is listed in an ordered form from high to low values 
and multiplicity of eigenvalues is permissible. 


Numerical Solution 


The solutions obtain from the system of equations (21) or (22) depend upon the 
selection of mean values associated with each energy interval. The selection of these mean 
values is determined by examing the numerical solution in certain special cases. We obtain 
a numerical solution of equation (8) in the special case given by g = g(E,x) — K Ee E ' T 
where K (# particles/cm 3 Mev) and T (Mev), are constants. We construct the solution 
over the spatial domain x > 0 and energy range .1 < E < 80 Mev. This domain is 
discretized by constructing a set of grid points Xi = iAx and Ej = jAE for some grid 
spacing defined by Ax and A E values. For i,j integers we define ujj = <p(xi,Ej), then 
the transport differential-integral equation (8) can be written in a discrete form. W e use 
the starting values uqj = 0 and vqj — 0. For the first step in Ax we write 

u i j = AxK Eje~ E]IT . (27) 


followed by the numerical calculation 



of 

( E')Te - r( - E '- E ^ 
1 _ e -(l-a)rE' 


u 




dE\ 


(28) 


when i = 1. After each numerical step the integrals of the type Vi j given by equation (29) 
are evaluated using Simpson’s 1/3 rule. We evaluate the equation (28) for j = 0, 1, . . we 
then use a two step algorithm in a repetitive fashion to advance the solution. For values 
of a near one the numerical solution of equation (8) requires that A E become small. For 


numerical accuracy we must have A.r << A E. The low energy spectrum then becomes 
difficult to calculate without special procedures, reference [1]. In this case we use a two 
step modified Euler predictor-corrector scheme, references [4], [5], which is defined by 

j? 

Second step: fij =v\ j + Eje 3 — cruij 

(uij + Axfij j = 0 

U2 ' J ~ \ i 1 + tiij+i) + Ax h tj j > 0 (29) 

Third step: f- 2 j —V 2 j + Ej€~ E: — cru- 2 j 

U 3 j -uij + 2 Ax f 2 j 
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The second step is an adaption of the Fredrichs method from reference [4]. The third step is 
a central difference second order step in Ax. After each 100 such applications of the above 
numerical two step algorithm, we apply the following stability correction as suggested in 
reference [5]. 

h.j =V3 ,j + Eje~ Ei -CT U3J 

1 (30) 

U3J =~ ( U3j + U2,j) + Ax / 3 ,j 


Recurrsive Solution 

In the special case g(E,x) = g(E), we assume a solution to equation (8) of the form 

oo 

<p(x, E) = ^ (j> n (E)f n {x) = (t>i{E)fi(x) + (j> 2 {E) f-2{x) H (31) 

n= 1 

We substitute this series into the equation (8) and obtain a solution by requiring that 


ME) 9 < fi/. /;w+»/iW=i (32) 

<t> n+1 (E)= f a (E,E')4>n(E')dE' n = 1,2,3,... f' n {x) + a/„(x) =/n-i(x) 

J E 

where the differential equations are subject to the initial condition that / n (0) — 0 f° r ^ n - 
Here the 4> n {E) terms are defined recursively and take a great deal of computational time 
for large values of n. The differential equations have the solutions given by the recurrsive 
relations 

/i(*) ) 


fn{x)= f In- \{u)e a{x u) du 

Jo 


(33) 


which are easily evaluated for as large an n as desired. We find numerically that \f n (x)\ 
decreases with increasing n when x is less than 1 and increases for x > 1 so that the series 
solution does not converge in this case. For |x| less than or equal to one we calculate the 
solution given by equation (33) for terms through n = 5 and n — 6 and compared them 
with the numerical solution. The mean values associated with the numerical solution of 
equations (21) and (22) where then adjusted in order that all solutions agree for this 
special circumstance. We then used these same mean values which where associated with 
numerical source terms as provided by the HZETRN code. 

Comparison of Multigroup and Other Solutions 

The numerical solutions and recurrsive solutions are compared with the multigroup 
solution for neutron penetration in lithium, aluminum and lead mediums for the case 
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n = 2. The results are illustrated in the figures 3,4, and 5. Excellent agreement is obtain 
in these cases. In these figures, the solid line represents the numerical solution. The 
circles represent the recursive solution and the triangles represent the multigroup solution. 
The various curves were calculated for depths of x = 0.1, 0.5, 1.0, 5.0, 10.0, 50.0 and 100.0 
centimeters. Similarly, the figures 6,7 and 8 illustrate the multigroup method in the case 
n = 10 for neutron penetration into lithium, aluminium and lead respectively. 

The multigroup method has a huge advantage in the computational time needed to 
calculate the solution. The multigroup method takes less than one minute of computational 
time while the other methods require many hours of computational time. 

Application for Al — H2O shield-target configuration. 

We now apply our previous development to an application of the multigroup method 
associated with an aluminium-water shield-target configuration. In particular, we consider 
the case where the source term g(E,x), in equation (8), represents evaporation neutrons 
produced per unit mass per Mev and is specified as a numerical array of values corre- 
sponding to various shield-target thicknesses and energies. The numerical array of values 
is produced by the radiation code HZEBIO, which is a modification of the radiation code 
HZETRN developed by Wilson, et. ah, reference [6]. The numerical array of values are 
actually given in the form g(Ei,Xj,yk) in units of particles /gm — Mev, where represents 
discrete values for various target thicknesses of water in gm/cm 2 , xj represents discrete 
values for various shield thicknesses of aluminium, also in units of gm/cm 2 , and E{ repre- 
sents discrete energy values in units of (Mev). We use these discrete source term values in 
the following way. We consider first the solution of equation (8) solved by the multigroup 
method with no target material i.e. all shield material with target thickness y = 0. We 
next consider the cases 2,3,... of discrete shield thickness X 2 ,x^,... and apply the multi- 
group method to the solution of equation (8) applied to all target material y > 0. For 
each aq-value considered, the initial conditions are obtained from the previous solutions 
generated where y = 0. This represents the application of the multigroup method to two 
different regions. Region 1 of all shield material and region 2 of all target material. We 
then continue to apply the multigroup method to region 2 for each discrete value of shield 
thickness, where the initial conditions on the start of the second region represents exit 
conditions from the shield region 1. This provides for continuity of the solutions for the 
fluence between the two regions. In this way we develop a series of graphs for fluence vs 
energy associated with various shield thickness. 
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In our application of the multigroup method we consider the case where the source 
term g(E , x), the scattering term f s (E, E') and cross section a(E) are all given numerically 
(obtained from the above HZETRN code). 

For a single shield material we solve 




pE/ a i 

4>(x,E)= / f Sl (E, E')<p(x, E f ) dE' + g(E,x) 

J E 

Integration of equation (34) from Ei to E {+ 1 produces 


(34) 




a(E)(f)(x, E) dE = 



f ai (E,E')4>(x,E')dE'dE + 



We define the quantities 


(35) 


<*>, = 



4>(x, E) dE 



(36) 


and interchange the order of integration of the double integral terms in equation (36). We 
then apply a mean value theorem to obtain the result 


d$i_ 

dx 


+ o 4q 



f Sl {E,E')dE <j>{x,E , )dE'+ 





lE^onE' 


f Sl (E, E') dE<t>(x , E') dE' + bi 


(37) 


where Ei < E' < Ei+\. The first double integral in equation (37) represents integration 
over the lower triangle illustrated in the figure 1(c). The second double integral in equation 
(38) represents integration over the upper triangle illustrated in the figure 1(c). Let 


f E' 

gi(E’)= / f Sl (E,E')dE 

JE=E t 


fEi+l 

and g 2 (E')= / f Sl (E,E')dE 

J E=aiE' 


(38) 


then employ another application of a mean value theorem for integrals to write equation 
(37) in the form 


d<|) 

— - + = gi{Ei + 6i(E i+ i - Ei))<t>i + g 2 {E i+ i + d-2{E i+ 2 ~ £i+i))$t-i + b i (39) 

dx 
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This produces the coefficients associated with the energy interval E{ to £;+ 1 given by 


an — gi — a and at,t + 1 = g 2 


( 40 ) 


needed for the numerical integration of equations (21). In this way the diagonal and off 
diagonal elements of the coefficient matrix in equation (21) are calculated. 

For a compound target material, comprised of material 1 and material 2, there are 
two values for a. A value a\ is selected for material 1 and a value c *2 is selected for the 
material 2 of the compound material. In this case the equation (34) takes on the form 


T + viE) 

OX 


<p(x, E) = 



fs t (E, E')<j>(x, E') dE'+ 



fsziE, E')4>(x, E') dE' + g(E,x) 


(41) 


where f Si and f S2 are scattering terms associated with the respective materials. These 
terms are calculated in the HZETRN code. We consider two cases. The first case requires 
that the E/a 2 line be above the E/a 1 line. (See figure 1(d)) The second case is where 
02 = 0 (the hydrogen case), and the limits of integration for the second integral go to 
infinity. We consider each case separately. 

For the first case we assume that ai > 02 > 0 and we select the exact energy spacing 
dictated by the E/a 2 line. We then proceed as we did using the single shield material. We 
integrate equation (41) from Ei to £i+i and interchange the order of integration on the 
double integral terms. Define bi = J E ' +l g{E,x) dE and obtain the equations 


d$i 

dx 


+ < 7 $i = /ll + 1 12 + 1 21 + I22 + 


(42) 


where now the /21 and I 22 integrals have, because of the exact spacings, the forms 

r&i+\ rE’ 

hi= / 


f S2 (E, E') dE(f>(x, E') dE' 


I Ei JE=Ei 

rEi+2 rEi+i 


(43) 


I‘22 


J E l 


i + l 


f E = Ol'2E > 


f S2 {E, E 1 ) dE<f>(x, E') dE' 


Defining the terms 



h2{i){E') 


Ei 


+ 1 


fsi (E, E 1 ) dE, i = 1,2 
fsi {E, E 1 ) dE i = 1,2 


E=a 2 E ! 
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and using the mean value theorem for integrals we obtain 

1 2 1 - /ii( 2 )(£; + Oi{Ei+i - Ei))$i and I 2 2 = ^2(2)(^i + ^ 2(^+2 - £*+i))*i+i 

where 0; and ^2 define intermediate energy values associated with the mean value theorem. 
The integrals In and I\ 2 are associated with integration limits (E,E/a 1 ) and energy 
intervals dictated by our selection of a 2 f° r determining the exact energy spacings. These 
integrals are associated with the trapezoidal area 1 and triangular area 2 illustrated in 
the figure 1(d). These areas are a fraction of the triangle area’s associated with the line 
E* — E/a 2 . These fractions are given by 


and we write 


j(Ej+i - Ej ) 2 - \(Ej + 1 - Ej/ai)(Ei+i - <*1^+1) 
fl Ei)(E i+2 -E i+ i) 

_ (E i+ i/ai ~ Ei+i)(E i+ i - Qiffj+i) 
h ~ {E i+ i - Ei){E i+ 2 - E i+l )) 


hi = f\h i(i) and In = 


The coefficients for our system of differential equations (22) are then given by 

ail =*1(2) + Mi(i) - ° ^ 

a 12 =h' 2 ( 2 ) + f 2 ^ 2 ( 1 ) ■ 

For the case 2, of hydrogen, a? = 0 and so one of the limits of integration becomes 
infinite. We let aq determine the energy spacing in this case. We again integrate equation 
(42) over and energy interval (E{, Fj+l) which is determined by the E = E /a i line. Lsing 
the definitions given by equations (36) we integrate the equations (41) over the interval 
(E{, Ei+i) and then interchange the order of integration in the resulting double integrals 

to obtain 


where 


—± + a<t>i = I* + I 2 + bi 

dx 


f E ' [ E f a .{E, E') dE<j>(x, E') dE 1 + [ f fsi(E, E 1 ) dE<p(x, E r ) dE' 
fEi J E—Ei JEi + i JE=ol\E' 


rEi + 1 pE' 


fEi JEi 


f 32 {E,E')dE<p(x,E')dE' + J2 


N rEi+j+i rEi 


f S2 (E,E')dE<p(x,E')dE' 


j = 1 d Ei+j 
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where for all N * greater than some integer N > 0 we know that f S2 (E,E ) will be zero. 
We define 

E 

h 3 (E') = f f Sl (E, E') dE , Ei<E' < E i+l 

JE i 
E 

h A {E') = I 1+1 f Sl ( E , E') dE, E i+ 1 < E' < E i+2 

JociE' 

r E‘ 

hs(E') = / f S2 ( E , E') dE, Ei< E' < E i+1 

JEi 

h§ (j) = / fs 2 (E, E r ) dE, Ei+j < Ej < Ei+j+i 
J Ei+j 

then we can write the coefficients associated with the system of differential equations as 


a iti =/i 3 + h$ - cr a i,i+3 —^6(3) 

&i,i+l =^4 + ^6(1) • 

a i,i+2 = ^6(2) a i,i+n ^6(n) 

In this way we generate a system of equations having the triangular form given by the 
equations (22). 

Various comparisons have been made to check the validity of the multigroup method. 
The figure 9 shows low energy neutron fluence vs depth for a shield-target aluminium 
water medium. Note the increase in low energy neutron production at the aluminium 
water interface at a depth of 100 5 /cm 2 . This is due to high energy neutrons colliding with 
hydrogen atoms. In these type of collisions the high energy neutrons give up over one-half 
of their energy, thus increasing the low energy neutron fluence. 

Using the source terms generated by the HZETRN program for an aluminium water 
configuration the figures 10 through 19 result for the evaporation neutron fluence as a 
function of energy for various shield thicknesses. The shield thicknesses in these figures 
are y = 0.0,0.3,1.0,5.0,10.0,20.0,30.0,50.0, and 100.0 g/cm 2 . The figures 20,21 and 22 
illustrate the comparison of the old HZETRN code results with and without the addition 
of the evaporation neutrons. These results are also compared with the Monte Carlo results 
for fluence associated with the February 1956 solar flare data. The multigroup method for 
the calculation of the low energy evaporation neutrons is computed much faster numerically 
than any of the previous approximation methods. The method also produced much more 
accurate results when compared to the Monte Carlo method, see foi example reference 
[157] and [158] of Appendix A. 
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Figure 2. Overlapping last position of multiple energy groups 
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Figure 4. Comparison of numerical solution with multigroup solution (n=2) for aluminum. 
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Figure 5. Comparison of numerical solution with multigroup solution (n=2) for lead. 
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Figure 6. Comparison of numerical solution with multigroup solution (n— 10) for lithium. 
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Neutron Fluence for Aluminum 



Figure 7. Comparison of numerical solution with multigroup solution (n=10) for aluminium. 
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Figure 8. Comparison of numerical solution with multigroup solution (n=10) for lead. 
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Figure 10. Multi-group low energy evaporation neutron fluence 

vs energy for target 0.0 gjcm 2 and various shield thicknesses. 
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Figure 11. Multi-group low energy evaporation neutron fluence 

vs energy for target 0.3 g/cm 2 and various shield thicknesses. 
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Figure 13. Multi-group low energy evaporation neutron fluence 

vs energy for target 3.0 g/cm' 2 and various shield thicknesses. 
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Figure 14. Multi-group low energy evaporation neutron fluence 

vs energy for target 5.0 gjcm 2 and various shield thicknesses. 
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Figure 15. Multi-group low energy evaporation neutron fluence 

vs energy for target 10.0 g/cm 2 and various shield thicknesses. 
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Figure 16. Multi-group low energy evaporation neutron fluence 

vs energy for target 20.0 g/cm 2 and various shield thicknesses. 
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Figure 17. Multi-group low energy evaporation neutron fluence 

vs energy for target 30.0 gjcm 2 and various shield thicknesses. 
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Figure 18. Multi-group low energy evaporation neutron fluence 

vs energy for target 50.0 g/cm 2 and various shield thicknesses. 
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Figure 20. Energy spectra of neutron fluence at 1.0 g/cm 2 depth of water slab 
exposed to February 1956 solar flare, as calculated by HZETRN 
with and without multi-group addition of low energy evaporation neutrons. 




NEUTRON FLUENCE FOR TARGET DEPTH 10 g/cm 2 OF WATER, NO SHIELD 





55 


Figure 21. Energy spectra of neutron fluence at 10.0 g/cm 2 depth of water slab 
exposed to February 1956 solar flare, as calculated by HZETRN 
with and without multi-group addition of low energy evaporation neutrons. 




NEUTRON FLUENCE FOR TARGET DEPTH 30 g/cm ! OF WATER, NO SHIELD 
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Figure 22. Energy spectra of neutron fluence at 30.0g/cm^ depth of water slab 
exposed to February 1956 solar flare, as calculated by HZETRN 
with and without multi-group addition of low energy evaporation neutrons. 














Appendix C 


Green’s Function Associated with the Transport of Light Ions in Matter. 
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Green’s Function Associated with the Transport of Light Ions in Matter. 


Basic Boltzmann equation 

^Change in ion flux ^ 
within a volume element 


/Gains within the\ 


V 


/ 


volume element 


V 


/Losses due to any\ 


/ 


nuclear collisions 


v 


J 


This gives the Boltzmann equation 




0j(x, fl, E) = Y l J dE ' f d ® T jkM*, SV, E' 

fa ^ j 


) ( 1 ) 


where 


(f)j(x, f2, E) 


E 

aj(E) 

Sj(E) 


Rj(E) 

j 

Cl 

rjk 


fl 

F 


X 

p 

X n 


is the flux of ions of type j moving in direction Q 

having units (#particles/cm 2 — sec — sr — Mev/amu) 
is the ion energy. (Mev/amu) 
is the atomic mass of the jth type ion (amu) 
is the macroscopic cross section, (cm -1 ) 
is the average energy loss per unit length or stopping power 
or linear energy transfer 4^. (Mev/cm). 
is the slowing down range for type j ions, (cm) Rj{E) = Jq 
is the ion type. 

is a unit vector in the direction of propagation. 

is production cross section of type j ions with energy E and direction Cl 
by collision with type k ions of energy E' and direction Cl' 
having units of (cm — sr — Mev/amu) 
is the outward directed unit normal to boundary, 
is vector to boundary point, (cm) 

is the position vector to arbitray point in region (cm) x = pCl + x n 

is the projection of x on Cl (cm) 

is the component of x perpendicular to 12 direction. 


The equation (1) is to be associated with the geometry of figure 1. 
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Figure 1. General Geometry for Boltzmann’s equation 
Multiply the equation (1) by Sj(E) and define the quantities 

4>j(x, fl, E) — S j(E)({>j(x, f2, E) 

Gj(xAe) = Sj(E)J2j dE' J dCl' r jk <Pk{x,Cl',E') 


to obtain 
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Aj dE V ’ 

9<Pi ■ 


(pj(x, Si, E) = G j(x, f2, E). 


( 2 ) 

(3) 

(4) 


Note that f 1 • V^- = is the directional derivative in the direction f2 and that 

dpj _ ddj dRj _ d(j>j Aj 
~QE ~ Tr~Te ~ dR~ Sj{E) 


so that the equation (4) can be written as 

$j(x,fi,E) = Gj(x,n,E). (5) 

Introduce the characteristic variables given by the transformation equations 


8 d 

b a(E) 

dp dRj v ' 


T)j = p-Rj(E) Z j = p + R J (E) 


( 6 ) 
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where p — f2 • x. Also introduce the variables 

Xj {Vj ■ £,j) — <Pj {p& + *?ni E) 

9j(Vj,Zj) = Gj{pCl + x n ,ft,E) (7) 

"AVjitj) = a j( E ) 




Figure 2. Geometry for characteristic variables 


By the chain rule we have 

Q<pj _ dtp d<pj 

dp drjj d£j 

so that the equation (5) simplifies to 
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and S=^(-D + ? j 


dRj drjj 


d Zj 


2^~ + <TjJ Xj(rij,Sj) = 9j{rij,Zj) 


(8) 


in terms of the new variables. This equation can be integrated using the integrating factor 

'1 W 


exp 
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to obtain 
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( 10 ) 
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where a is any real number. Consequently, the solution to the equation (4) can be written 


as 
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J f(v'\ p ~ Rj (E))Gj + V / )n + gn,n,fl ~ 1 ( ~" 2 )J dif 


(ii) 


where 


/(a, p - Rj(E)) = exp 
From equations (6) we find that 


1 rp-Rj(E) c . _ n / 

2 / 


( 12 ) 


2 P = Vj + 0 and 2Rj{E) = - pj 


(13) 


so that when r\' — a we will have p = A(a + Observe from figure 2 that along the line 
of integration we will have = constant. The value of a is selected such that pCl +x n = f 
is a point on the boundary. Thus, the vector j Cl + x n = f dotted with C, gives the 


value 


a = 2Cl-f -ij = 2d- p- Rj{E) (14) 

where d - ff ■ F. Note that when E — E 1 and pj = p' we have from equation (13) that 


2 R j {E')=i j -r! 


(15) 


or 
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dp' (16) 


and similarly by changing symbols when 
E" = R 


" - "7 1 ( P + we have d£" 
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We examine the limits of integration in equation (11) and observe that when p 1 
have 

2 Rj{E') = p + Rj{E) - a 

and from equation (14) we have 


(17) 


a we 


(18) 


2d — p + Rj(E) + cl. 


(19) 
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Adding the equations (18) and (19) we find 


Rj{E')+d = p + Rj{E) 


( 20 ) 


or 


E' = R~ l (p — d + R-j(E)). 


( 21 ) 


Next we examine the lower limit of integration and find that when rj 1 = p — Rj(E), then 
2 Rj{E') = p + Rj(E) — p + Rj(E) implies that E' = E. In the second term of equation 
(11) when r /' = a we again find that E" = RJ l (p + Rj(E) — d) and when rj" = p — Rj(E) 
then E" — E . Also, 

+ v") = + ti ~ 2 Rj(E")) =Sj- Rj(E") = p + Rj(E) - Rj{E"). 

Consequently, the equation (11) can be written in the form 


(x, n, E) = Fi (E, RJ 1 (R 3 (E) - d + P)4 >j (f , n, RJ 1 (Rj ( E)+p - d)) 

rR- 1 (R j (E)+p-d) _ 

J El (E, E")Gj ((p + Rj ( E ) - Rj (E"))Q + x n , fi, £") ^§77} dE r 


+ 


where 


Fi(El, E 2) = exp 


_ r^n dE ' 

J El Sj(E') 


Define the nuclear survival probability (reference Wilson 1977) as 


Pj(E) = exp 


-f 


Iq^- dE' 

Sj(E') 


then the equation (23) can be written as 

Fi(E u E 2 ) 


PjjE 2 ) 


( 22 ) 


(23) 


(24) 


Then from equation (22) we can write the solution to equation (3) in the form 

+ S/ ' dE> ' s^(E)P^(E) / dE " / E")<t>h (x + (R 7 (E) - Rj(E'))Q, ft' , E") 


( 25 ) 


where Ej = Rj l (p + Rj{E) — d), x = x n + pCl and E' and E" have been interchanged. 
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In the one-dimensional straight ahead approximation Cl is a unit vector in the direction 
of x with p = x, x n = 0, r)j = x - Rj(E),£j = x + Rj{E) and f = 0. (i.e. the origin 0 
moves to the boundary x = 0). The equation (25) then reduces to 


, (r n _ S j {E j )P j (E j ) 


+ E f ’ dE 'CrtCinh I iE"<, jk (E’,E")0 k (i + Rj(E) - Rj(E'), E") 


(26) 


where Ej is determined from x and E such that 


Ej = RJ l {x + Rj{E)). 


(27) 


The solution given by equation (26) 


can be expressed in terms of Green’s function as 



Gjk(x, E , £o)0fc(O, Eq) dEo 


(28) 


where <f>k(Q,Eo) = fk{E o) are boundary conditions. Substituting the assumed solution 
given by equation (28) into equation (26) we obtain 


w_ » r°° 

Y / Gjt{x,E t 

t 

*?/. 


Eo)4>e (0, Eo)dEo = U 




Sj(Ej)Pj{Ej) 

Sj(E)Pj{E) 


Gjt( 0, Ej, Eq)4>( (0, Eo)dEo 



G ke (x+R ] (E)-R j (E'), E", E 0 )4> t (0, E 0 ) dE 0 . 


Note that when l = m we can equate like coefficients and find that Gj Tn (x,E,E o) must 
satisfy the integral equation 


I T _ Sj(Ej)Pj(Ej) ^ m t> -o \ 

Gj m (x, E, Eq) — ^rT^rR~rE\ Gj m {0, Ej , Eq) 


+ 


Sj (E)Pj(E) 


k>j 


r ngTTtt f dE” <T jk (E',E")G km (x + R 3 (E) - Rj(E'),E",E 0 ) 


(29) 


subject to the boundary condition Gj m ( 0, E, Eo) = 6jm6(E — Eo), where the value for Ej 
is determined from the inverse relation Ej = RJ l {x + Rj(E)). The Gjm terms are written 
using the Neumann expansion as a perturbation series 

CO 

G jm (x,E,E 0 ) = Y G jm( x ' E ' E °) ( 30 ) 

i=0 
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with leading term 


£.(0)/ r p p R j( E j) R j( E j) c r/p p \ 

^ jm( x > £o) — S (E)P (E) Eo). 


(31) 


with Ej = Rj 1 (i + Rj(E)). Note that when x = 0 we have Ej — E so that G^(0, E, Eo) 
satisfies the above boundary condition. The higher order terms are determined from the 
recursive definition 

G^ l \x,E,Eo) = 

j p> AjPj{E ) f ipff / pt „ /pi p /p/-, p// r i (*^) 

E dt 'S ~{E)P j {E) J E , ^ j>C ^ E >E > G km( X + R j( E ) ~ R j( E )’ E > E 0)- 


?/. 


and must satisfy the boundary conditions (x ) E, Eo) — 0 for n = 0, 1, 2, . . . . In the 

special case n — 0 the equation (32) reduces to 

GV(x,E,E 0 ) = 

V f Ei dF'-MA £ L f^zr" /p/ F //x SkjE^PkjE^ _ (33) 

\ Je S j (E)P j (E) J e 1 ’ h k {E")P k (E») Skm ^ k 0) 

where R k (E' k ) = x + Rj{E) - Rj(E') + R k {E"). (i.e. treat x + Rj(E) - Rj(E") as an 
x* value. See for example equation (27) .) Again we observe that when x = 0 we have 
Ej = E and so the boundary condition at x = 0 is satisfied. 

Cross Section assumption 1 

For interactions dominated by perpherial processes we use 


a jm (E',E") = a jTn (E")6(E' -E”) 
so that the equation (33) becomes 

Gfl{x,E,E 0 ) = 


(34) 


k 

where 


J F « -j-jy-, i’°°,p// T (P»\C(P> pH\ S k{E , k )P k {E , k ) , (3o) 

P-Je dE —E)rAE)J B , iE )S{E ~ E 


E'k = + R,(E) - Rj(E') + R k (E'j). 


(36) 


We integrate with respect to E" and observe that the only nonzero term occurs when 
E" = E' . This gives 


ntl) (r p p y = V I** dF' ±j£ilEl - rr / ^fc( £ fc) p fc( g fc) c rr / _ F y s , 

Jm ' ’ ’ V Je E Sj{E)Pj(E) ^ S fc (F / )Pfc(F / ) ^ * 


(37) 
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where 


Rk(E' k ) =x + Rj(E) - Rj(E') + Rk(E') 
Rj(E') - R k (E') =x + Rj(E) - R k (E' k ). 


(38) 


We know that i /jRj(E') = iy k R k (E / ) so that the above can be written as 




- 1 R k {E') -x + Rj(E) - R k (E' k ) 


or M - -ij RjiE') -x + Rj{E) - R k (E' k ) 


Thus, we can write 


Rj(E') = ] ~~~ " ' | (x + Rj(E ) - R k {E k )) 

K - "jl 

or R k (E') = - Vi (x + Rj{E) - R k (E' k )). 

K ~ Vj\ 

Differentiate the equation (39) with respect to E k to obtain 

R' k (E')dE' = - - Vi - (-R' k (E' k ))dE' k or dE' - SkiE ^^> 

\Vk Vj I 

The equation (35) can then be written as 

"El 


\ Uk Sic (E 1 ) 


e, *) = £ / “ 

k ^fcl 


Pk(E') 


(39) 


i±dE' k (40) 


Pk{E 'b± 6 {E' k - E 0 )8 km ( 41 ) 


The only nonzero contribution comes when k = m and E k — Eq and so equation (41) 
reduces to 


/■ ]\ ( hjm E , Eq , £/ ) if (i?m (So) Rj (E) <C. Rm {E q) % 

J (42) 

1 0 otherwise 


where 


and 


h jm (x,E,E 0 .E') = JkElEl g (E') ^ Pm(Eo) 

J V ' Sj{E)Pj(E) Jm{ ’ \v m ~ vj\ Pm{E') 


E' = R~ l 1 Um 
J ' l^m 


7 [* + /2 i (E) -fl m (£o)]J • 

-^1 / 


(43) 


(44) 
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That is, when E' k = Eq and k — m we have from the equation (3b) that 

R m {E 0 ) =x + Rj{E) — Rj{E ) + Rm(E ) 

Rj{E') - Rm{E f ) —x + Rj{E) - R m {E 0 ) 

Rj(E') =x + Rj{E) - Rm{Eo) 

Rj(E') = , Um — r(x + Rj(E) - Rm(Eo)). 

Wm-Vj I 

Also from the transformation equations (13) the r]k, ^kiVj^j variables are related thiough 
the range scale factors Vj and Vk, where VjRj = VkRk • This produces the relations 


Vk~ Zk = -2Rk = -2 -i-Rj - ~ tj)- 

Vk Vk 


Then from the equations 


tj + Vj + V k 


Vj ~ tj - — {Vk ~ tk) 

v j 

we find that by adding the equation (45) and (46) we obtain 

*lj = ( 

and subtracting (46) from (45) we obtain 




Vj. 


2*; 


Vk 


Vj 


Interchanging j and k in the equations (47) and (4b) we find that 

. V ✓ \ 

Vj + Vk 


Vk 

Zk 


2v k , 
Vk - Vj 

‘2vk 




V, + 


2v k . 
Vk + Vj 

2vk 


ti- 


HS) 

(46) 


Vk 

1 - — 


(47) 

v ij 

/ 


Vk ^ 
1 + — 

) ik- 

(48) 

V j) 

r 



Then when 77 j is a value r/ lying between the constants — £j and +£j , (See Figure 2(b)), 
we will have 

/ . -L \ 

V + 


Vk = 


Zk = 


Vj + v k 


2vk 

Vk - Vj 

2 vk 


V + 


v _ ’Jc v j 

2 v k , 
Vk + Vj 
2 vk 


£? 

tj- 
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Changing k to m we find 


£m 


I'm ~ Vj 
•2 I'm 


V + 


( "rn + Vj \ 

V 2z/ m 


Note that the boundary condition Gj m (0, E, Eq) = 6j m 6(E — £o) can be written in the 
form 


C jm (0, E, E 0 ) = 6 jm 6(R- l (Zj) - E 0 ) = S jm 6(^ - Rj(E 0 )) = <5(£ m - fl m (E 0 )) 
so that when = R m (E 0 ) we have 




2 I'm 


^m(^o) 



(49) 


Using the equations (6) and (49) we now calculate the inequality which occurs in the 
equation (42) . From the equation (10) , with a = — £, we have the inequality — < rj' < rjj 

which implies 


-x Rj(E) < 


2z/ r 


<r)' <r) j 
Vm E V 


Vrn Vi 


■R m (Eo)— 


—x < 


2i/ r 


I'm - V 


■Rm{E 0 ) + Rj(E)~ 


, u rn - Vj 
I'm + V ' 


2. (x + Rj(E)) < X - i?j(E) 


I'm ~ I' 


2- (x + Rj{E)) < X 


I'm + I'? 


X — X < 


‘2v m R Tn (Eq) 


'3/ I'm ~ I'j 

I'jX <! v rn R m (Eo) vjRj(E) < v m x 


+ 1 - 


I'm + I' 


^771 


- U;(£)<x + 


J-'m “t" 


V 7 


Vjjx R m ( £*q ) > VjnX 

VmRm(Eo) - VjX > UjRj(E) > v m R m (Eo) — V m x 

—R m (E 0 ) - x >Rj(E) > — (R m (E 0 ) - x) 
v J v i 

— (Rm(Eo) - x) <Rj(E) < —Rm(Eo) - X 
U 3 I 'j 


Cross Section assumption 2 

We start with equation (33) and assume aj m (E',E") has a Gaussian distribution of 
the form 


G jm (E',E") = a jm (E") 


j\jm \/27 


exp 


(£' ~ E" Z W 
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we can then write the equation (33) in the form 


where 


pEj poo 

GyJ l (x > E,E 0 )= dE' dE"F jk 8 kr 

JE JE' 


_ AjPpEy , SjyEjyPkjEy) , _ 

3 Sj(E) Pj(E) Jk{ ’ ) S k (E”)P k {E") { k 0) 


(50) 


(51) 


with 


E' k = R-^x + Rj(E) - Rj(E') + R k (E")) (52) 

The integration of (50) is over the region illustrated in the figure 3 in the limit as T — ► oo. 
In expanded form the equation (50) has the form 


where 


G { yl{x,E,E 0 ) = 




gfc(flfc)iMsfc) 

S k {E")P k {E") 


t>km8(E k ~ Eq). 


®jk ) 


Ajks/m 


exp 


{E'-E"-e jk )< 


( 53 ) 



Figure 3. Limits of integration for Green’s function term. 
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Note that for the first integration in the E" direction we have E' is a constant. 
Consequently, we let 


r = 


E"+e jk - E' 


with dr = 


dE r 


\/2A jk \/2A j k 

The equation (53) can then be written in the form 
G^(x,E,E o ) = 


(54) 


r E i r°° 

L 

V2A jk 


AjPjjin _ r2 5 fc (gj()P fc (g^) c , (55) 

d yf* S j {E)P J (E ) Jk( ) S k (r)Pk(r) ^ ( fc ' 


where r = \/2Aj k r - ej k + E' and 


E' k = R~\x + Rj(E) - Rj(E') - R k (f)). 


(56) 


This integral can be simplified by using one of the mean value theorems for integrals and 
written as 


G^(x,E,E 0 ) = 


L 


AjPjiE 1 ) _ S k {E' kt )P k (E' kt ) 


dE' J —^ — —a, k (f .) - 

E 2S 3 {E)P ] {E) > kK 


S k {f.)P k {f.) 


S km S(E' kt - E 0 ) 


_ 2 _ r°° 

y/* J-ALh— 


e~ r dr. 


(57) 


with r* = \/2A j k r* — ej k + E' and 


E'k. = ** (* + Rj(E) - Rj(E') - Rk(r ,)) 


(58) 


where r* is some mean value in the interval ( -^A— , oo) and when E' k * = Eq, then E' is a 
solution of the nonlinear equation 


R k (E 0 ) = x + Rj{E) - Rj{E') - R k (V2A jk r* - e jk + E') 

provided E < E f < Ej. Consequently, we can write 
gM(x,E,E o) = 

1 AjPj(E f ) f \ Sm{Eo)Pm( Eq) p r ( *Jrn \ ]f r / p/ y r. 

I 5 ; mrr l E) a :rn (r. ) 3 ^ 7 y p m (P . ) ertc J it E < E < E } 

0 otherwise 


(59) 


(60) 


where E 1 is a solution of the nonlinear equation (59) , r* is some mean value and erfc is 
the complimentary error function. 
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Figure 4. New limits of integration for Green’s function term. 

Another viewpoint 

By interchanging the order of integration in equation (33) we obtain the limits of 
integration illustrated in the figure 4 and the equation (33) can be written as 

G^{x,E,E 0 ) = J 3 dE" J dE'F jk 6 km + dE " J e dE'F jk 6 krn . (61) 

Observe that along the line £’ // =constant, we have from equation (52) that 


dRk(E' k ) dE ' k _ 
dE' k dE ' 

Ak dK _ 

S k {E' k ) dE' 
or dE ’ = 


dRj{E') 

dE' 

Aj 

Sj(E') 

AmSj(E') jr , 
AjS k (E' k ) fc ‘ 


Hence, when k = m and 6 km = 1, the equation (61) reduces to 





A m Sj(E') 


dE' + 


lim 
T — oo 





AmSj(E') 

AjS m (E' m ) 


dE' 


( 62 ) 
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The limits in the above equation are determined as follows. Observe that when E E 
and k — m the equation (52) gives 

Rm{E' m ) =x + Rj(E ) - Rj{E) + Rm(E") (63) 

or E' m =R^(x + Rm{E”)) = E m i 
and when E' = Ej and k = m the equation (52) gives 

R m {E' m ) = x+ Rj(E) - Rj(Ej) + R m (E") 


. But Rj{Ej) = x + Rj{E) so that E' m = E" . Also when E ! = E" and k - m we obtain 
from the equation (52) that 


R m (E'J =x + Rj(E) - Rj(E ") + R m (E ) 

/M i/jRj(E") 

R m (E' m ) =x + Rj(E) - Rj{E”) + 

v m 


Rm(E’ m ) =x + Rj(E) + - l) Rj(E") 

or E’ m =R~H* + R,{E) + - i) Ri(E")) = E m 3 


with 

E' = RJ l (x + Rj{E) - Rm{E' m ) + Rm(E"))- 

Using the properties of the Dirac delta function we find that the only nonzero contribution 
to the integral dE^ occurs when E' m = Eq. In this case the integral gi\en by equation 
(62) simplifies to 


Gfl(x, E, Bo) 


jm ' 

f B ,„„A m S 1 (E')P,(E>)_ , fin(Eo) f 

J Ei dE SJSW^E T Jml ' S m (E")P m (E") 11 

[ T „A m S,(E’)P,[E') , Pm(Eo) 

+ & J E dE S J (E)P J (E)~ aim{E ' ] S m (E")P m (E") } 


(64) 


where 

E' = Rj\x + Rj{E) - R m (E 0 ) + R m {E ")) 


(65) 


and 



f i 

if Em3 < E 0 < Emi 

/ 1 = 

(o 

otherwise 


r i 

if E" < Eq < E m i 

f‘2 = 

l o 

otherwise 


( 66 ) 

(67) 
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